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^SJ ■ Abstract 



We perforin Monte Carlo simulations using the Wolff cluster algorithm of the XY model on both 
fixed and dynamical phi-cubed graphs (i.e. without and with coupling to two-dimensional quantum 
gravity). We compare the numerical results with the theoretical expectation that the phase transi- 
tion remains of KT type when the XY model is coupled to gravity. We also examine whether the 

/r\ " universality we discovered in our earlier work on various Potts models with the same value of the 

r^ , central charge, c, carries over to the XY model, which has c = 1. 
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1 Introduction 

In the past few years the work of KPZ in the hghtcone gauge LiouviUe theory ||l[ and DDK in the 
conformal gauge g| has aUowed the calculation of the critical exponents for c < 1 conformal models 
coupled to 2d quantum gravity. These conformal field theory results are backed up by the exact solution 
of the Ising model, with c = 1/2, using matrix model methods in p], which gives identical exponents to 
those calculated in 0] , ■ Numerical simulations of the Ising model on both dynamical triangulations 
and dynamical phi-cubed graphs [El, which we would expect to discretize the continuum theory coupled 
to gravity, also give good agreement with these analytical predictions. We have recently simulated the 
Q = 3, 4 state Potts models on dynamical phi-cubed graphs, with c — A/b, 1 respectively, and again found 
good agreement with the conformal field theory predictions B (in these cases there are, as yet, no exact 
solutions available). 

All of the analytical approaches to 2d quantum gravity, whether by continuum LiouviUe theory or 
matrix models, break down at c = 1. We can see this quite clearly in the KPZ/DDK results for the 
"dressed" conformal weight A of an operator with weight Aq before coupling to gravity 



A-Ao = -yA(A-l), (1) 



where a is 



a = - — {^/W^-^/T^. (2) 

We thus get nonsensical complex weights for c > 1. Similarly, d = multimatrix models and the 
multicritical points of the d = Q single matrix model appear to give access only to c < 1. There have 
been various speculations regarding the nature of the "barrier" at c = 1 M and suggestions on extending 
the validity of the LiouviUe theory to c > 1 H , although it is interesting to note that both numerical B 
and analytical p] work with multiple Potts models (with c > 1, at least naively) have failed to see any 
clear signal of pathological behavior. 

In spite of this caveat it is clear that c = 1 represents a particularly interesting case. In the matrix 
model approach we can construct a c = 1 model by considering matrices in d = 1 which gives the partition 
function 

(3) 



Z ^ I L>$exp -13 I dtTr ( -4^ + [/($) 



where $ is an A'' x A'' hermitian matrix. This model was solved in the N ^ oo limit in pU and shown to 
be equivalent to c = 1 matter coupled to 2d quantum gravity in |ll|. More recently the double scaling 
limit of the model was taken in [O . As the model's continuum limit is an uncompactified boson it is 
natural to try and extend the analysis to the case of a compactified boson, which gives us the XY model 
coupled to 2d quantum gravity. 

The standard XY model on a regular lattice (i.e. not coupled to 2d gravity) with a partition function 
of the form 



z = l[ fde,expipJ2' 



cos(0, - e^) (4) 



displays an infinite order "Kosterlitz-Thoulcss" (KT) phase transition caused by the liberation of vortices 
[fi3[ . At this transition the correlation length diverges as 



and the susceptibility as 



X = .4, »p (p.1^) . (6) 

where the critical exponent ly is predicted to be 1/2 p3|. T is the temperature defined with A; = 1 so 
that P = 1/T. This KT theory also predicts the correlation function critical exponent rj = 1/4, where rj 
is given by 

x = ce-"- (7) 



In spite of some earlier controversy |l4| , apparently due to difficulties in distinguishing KT fits from 
second order ones, more recent numerical simulations support these results [|l5| . 

Returning to the XY model coupled to 2d gravity, the compactified model was solved in jlj] for 
a singlet sector in which the angular integrations in the matrix integral were unimportant. This was 
identified with the vortex free sector of the XY model. The adjoint sector was then incorporated into 



the model in ||l7| and identified with a single vortex/anti-vortex pair. The predictions of |I7 for the XY 
model coupled to 2d quantum gravity were, in direct analogy with the standard XY model, a KT phase 
transition at (3c — Stt. It is interesting to note that the phase transition is predicted to remain the same 
when coupled to gravity, unlike the case of the Ising and Q = 3, 4 state Potts models where second order 
phase transitions are weakened to third order. In this paper we simulate the XY model on dynamical 
phi-cubed graphs 

^^ = E n / '^^^ exp I -/3 J2 ^if^ cos(0, - ^,) I (8) 

where G^ is the connectivity matrix of the graph, in order to see how these predictions tally with 
"experiment" . In B we found universal graph properties for various Potts models, depending only on c, 
so we also measure these here to see how the XY model compares with the other c = 1 models we have 
investigated. 

2 Fixed Simulations 

As a check of our program we first simulated the XY model on a fixed random lattice coming from a pure 
two-dimensional quantum gravity simulation p8[ . These lattices are phi-cubed graphs with the topology 
of a sphere, so each point has exactly three neighbors. It is currently unclear whether it is possible to 
take the continuum limit on a single graph of this form |19|] because of its highly irregular nature, but 
we are principally interested here in using the simulation as a test to compare with our simulation on 
a dynamical graph so we do not concern ourselves with this problem. In any case, a simulation of the 
Ising model on such a fixed graph gives results that are compatible with the standard square lattice 
Onsager exponents |2y] . Our simulations use the Wolff cluster algorithm [|TJ to update the spins which is 
much more effective in ameliorating critical slowing down than traditional (local) Monte Carlo methods. 
XY models on random triangulations have been simulated in the past using the Langevin method |p2| 
and the Monte Carlo Renormalization Group method |23[| . The random lattices used in these previous 
simulations were obtained by linking up into triangles points distributed at random in a unit square 
with periodic boundary conditions, which gives a much less disordered graph than that coming from the 
two-dimensional quantum gravity. 

As this is merely a test simulation we used only three small random graphs with N = 100, 500 and 
1000 points. At 27 values of (3 between 0.1 and 3.0 we ran 10000 sets of cluster updates to equilibrate 
the XY spins and then measured the energy E, specific heat C, susceptibility x ^-i^d correlation length 
^ over the next 50000 updates. (The magnetization AI is a vector for the XY model with components 
(^jCOs6'i,^jSin0i) so one usually measures instead the susceptibility which is defined as its square 
X =< M ■ M > /N. The energy, specific heat and correlation length are defined in the usual way.) 
Depending on f3 each set of cluster updates consisted of up to 100 "hits" of the Wolff algorithm, this is 
necessary because as /3 decreases (the temperature increases) the size of the cluster built by the Wolff 
algorithm decreases therefore more applications of it are required to update the same number of spins. 
To verify that the Wolff algorithm had overcome critical slowing down we measured the autocorrelation 
time r for the energy, and found that it increases slightly as /3 increases but does not increase with N, so 
critical slowing down has essentially been eliminated. The autocorrelation time in the data is always less 
than 10 updates and we bin this data in bins of size 1000 for statistical analysis so our error bars can be 
trusted as coming from uncorrelated samples. 

In our earlier simulations of single and multiple Potts models coupled to 2d quantum gravity |q, p| , we 
were able to make use of Binder's cumulant to determine the critical inverse temperature (3c and critical 
exponent v separately. Unfortunately this method cannot be applied to the XY model so we must obtain 
these from the four-parameter fits in eqs. ra,pl As is well known this is rather difficult - even with very 
high statistics data on large lattices - hence the controversy mentioned earlier |l4, O . Things are made 



worse by the behaviors of the correlation length and susceptibility in the XY model: they both diverge 
according to eqs. ||J| for T > Tc {(3 < /3c) but are infinite (on an infinite system) for T < Tc {(3 > pc)- 
Therefore there are no peaks at Tc and on a finite lattice ^ and x just keep increasing as T decreases, 
with any jump at Tc rounded away by finite-size effects. The specific heat does have a peak, which is 
independent of system size. Fig. 1, but this occurs at a lower value of (3 than f3c so is of no help in 
determining Pc- Hence from our data it is to all intents and purposes impossible to determine both Tc 
(Pc) and z/. The only way to proceed is to assume that v — f /2 as predicted by KT and attempt to fit 
to Tc using the resulting three-parameter fits: 



?-^cexp^^^=^j, X = A,exp^-^^^J. (9) 

To get an idea of what Pc (and thus Tc) is approximately we do some finite-size scaling using eqj^. In 
finite-size scaling relations one usually has some quantity scaling as a power of the linear dimension, L, 
of the d-dimensional system with N — L points. Here we do not know d, the internal dimensionality of 
our system, a priori so must write instead Nd . The method we shall adopt was used in [M to estimate 
the exponent rj knowing Pc', here we shall use it the other way around, assuming that 77 takes its KT 
predicted value of 1/4. This is a reasonable thing to do because the XY model has a line of fixed points 
for P > Pc along which x and £, diverge according to eq.|7|, however, rj depends on P and in fact the KT 
analysis predicts rj ~ 1//3, as p ^ 00. Therefore we only obtain the value ry = 1/4 at Pc- The method 
works as follows. For P > Pc the correlation length diverges and therefore becomes the system size Nd 
on our finite lattices. Thus eq.0 becomes 

X - CiV(2-")/''. (10) 

Therefore if we plot [ln{x) — ln{C)]/ln{N)] versus P then we should obtain (2 — ri)/d for P > Pc- At Pc 
where rj — 1/4 and we take d = 2.6 from below, the expected value is 0.673. We tune the value for ln{C) 
until all the points for P > Pc from the different N fall on a universal curve. The resulting plot (with 
ln{C) = 1.0) is shown in Fig. 2, from it we estimate Pc = 2.3(1). 

We now bite the bullet and attempt to fit to ^ and x- For ^ fitting the data between T — 0.5 and 
0.9 yields Tc = 0.36(1) with a x^ /dof = 7.9. At T = 0.5 the correlation length has grown to 10 however 
which is rather large for the N — 1000 system. Therefore we restrict the fit to the 7 points between 
T = 0.6 (where ^ « 6) and 0.9 to obtain Tc = 0.39(2) with a x^/do/ = 0.3. Taking the difference as a 
measure of the error we assume Tc = 0.39(3) so that Pc — 2.6(2). This value is rather high compared 
with the 2.3 we obtained from finite-size scaling but, from dynamical results below, we know that this is 
a result of having too small a lattice - on a larger lattice it would come down to at least (most) 2.4. We 
cannot get nearly as good fits from x, fitting the same two temperature ranges yields Tc = 0.37(2) and 
0.43(2) with enormous x^ /dof^s. This fortuitously gives Pc — 2.3(3) but again for larger lattice sizes this 
value will decrease. Even for the large-scale simulations of the regular XY model ||l^ fits to x are almost 
always worse than fits to ^. From this the best we can conclude is that Pc = 2.3(1) for the XY model on 
a fixed random lattice coming from quantum gravity. 

3 Dynamical Simulations 

For the dynamical simulations, we started from a random phi-cubed graph coming from pure two- 
dimensional quantum gravity but updated it during the simulation using the standard "flip" move with 
the detailed balance condition that the rings at either ends of the link being flipped have no links in com- 
mon. This check eliminates all graphs containing tadpoles or self-energies. Therefore the Monte Carlo 
update consists of two parts: firstly the XY spins are updated by a set of Wolff "hits" then Njnp links of 
the graph are picked randomly and flipped. After testing various values of NfHp to ensure that there were 
enough flips to make the graph dynamical on the time scale of the XY spin updates we set NfHp — N for 
all the simulations. We ran on iV = 100, 500, 1000, 2000 and 5000 at 27 values of P between 0.1 and 3.0, 
doing 10000 updates for equilibration and 50000 for measurement as for the fixed case. In addition to 
the standard thermodynamic quantities for the spin model: E, C, x and ^ we measured some properties 
of the graph: acceptance rates for flips, distribution of ring lengths and internal fractal dimension d. We 
must measure d because we do not know a priori the internal dimensionality of our system. In fact. 



numerical simulations for pure 2d quantum gravity |25 and for 2d quantum gravity coupled to Potts 
models ^, g| yield d « 2.7; and an analytical calculation predicts that d lies between 2 and 3 for pure 2d 
quantum gravity p6[ . 

We measured the autocorrelation time r for the energy for N — 1000, 2000 and 5000 and found that 
it increases slightly as (3 increases and with N near the phase transition, though never exceeding 9. At 
the critical point /3c, where the correlation length becomes the size of the system, we can fit 



N^. 



(11) 



to extract an estimate for the dynamical critical exponent z. If we take our best estimate of Pc from 
below and assume that d — 2.6 then we obtain z = 0.8(1). This is larger than what we found for Potts 
models g, g], but still a lot less than traditional (local) Monte Carlo methods which have z « 2. 

As for the fixed simulations we have the same problems in fitting Tc from the correlation length and 
susceptibility. We restrict both fits to 8 points between T — 0.55 and 0.9, so that ^ < 6, and obtain the 
values for Tc listed in columns 2 and 3 of Table 1. We see that, for ^ at least, the Tc values increase 
slightly as N increases. Taking the difference between the results on iV = 2000 and N = 5000 as an 
indication of the error, we thus obtain /3c = 2.4(2) from ^ and /3c = 2.0(2) from x- The fits for ^ yield 
X^/do/'s of around 1, whereas those for x again are huge. Thus we take 2.4(2) as the best value for /3c. 
This value is confirmed by our finite-size scaling plot of [ln{x) — ln{C)]/ln{N)] versus /3 in Fig. 3 which 
yields /3c = 2.4(1). Finding a /3c for the dynamical case that is larger than the fixed case satisfies naive 
expectations, since this means that it is more difficult to "freeze" the dynamical system (so we must go 
to a lower temperature). 



N 


Tc from ^ 


Tc from X 


d 


1000 


0.38(1) 


0.50(6) 


2.13(2) 


2000 


0.39(1) 


0.51(6) 


2.25(2) 


5000 


0.41(1) 


0.51(5) 


2.42(1) 



Table 1: Fitted values of Tr and d. 



Again the specific heat does have a peak occurring at a lower value of /3 than Pc- However, it now 
appears that the peak is growing with system size (which would be the behavior expected from a second 
order phase transition), as we can see in Fig. 4. We think that this is an artifact of not running long 
enough to obtain enough bins in order to estimate the specific heat accurately from the fluctuation in 
the energy: 

P' 



C=^{<E'>~<E>'). 



(12) 



Evidence for this is that by using the other way to calculate C ~ numerically differentiating the energy 
- we obtain curves, which we also plot in Fig. 4, that look exactly like those from the fixed simulations, 
where both methods for obtaining C agree. 

Turning now to the graph properties. We measure the acceptance rate for the flip move to confirm 
that our graphs are really dynamical. The flip can be forbidden either from the graph constraints coming 
from the detailed balance condition or from the energy change of the spin model, so we can decompose 
the flip acceptance rate into two parts: AL - the fraction of randomly selected links which can be flipped 
satisfying the graph constraints; and AF - the fraction of links satisfying the graph constraints which are 
actually flipped, i.e. pass the Metropolis test using the XY model energy change. AF and AL are shown 
for N — 2000 in Fig. 5. We see that they both have dips at some P < Pc- 

The distribution of ring lengths in the graph is the discrete equivalent of the distribution of local 
Gaussian curvatures in the continuum. The minimum possible ring length is 3. If we also plot the 
probability (fraction) of rings of length three (PR3) in Fig. 5 to compare with AF and AL we see that 
it has a peak very close to the dip in AL. This is reasonable since both PR3 and AL depend only on the 
graph, whereas AF depends on the spin model. Interestingly the peak in PR3 and dip in AL occur close 
to the peak in the specific heat, at a value of /3 lower than /3c. In Fig. 6 we compare PR3 with what we 
obtained for single Potts models for A^ = 2000 in order to see if the universality we discovered carries 
over to the XY model (note that we have had to scale /3 for the XY model by a factor of two to correct 



for differing definitions in the actions) . f f it did we would expect PR3 for the XY model to look like that 
for the Q = 4 Potts model as they both have central charge f . ft is rather obvious that PR3 for the XY 
model is significantly different, in fact its maximum value of 0.2f 85 is rather close to that of the fsing 
((5 = 2 Potts) model (0.2179). However its value at the phase transition 0.2160 lies between those for 
the Q = 2 and (5 = 3 Potts models (0.2151 and 0.2167 respectively), and a long way from that for the 
(5 = 4 model (0.2190), so must conclude that there is no universality between the Potts and XY models 
coupled to quantum gravity. 

Finally we measured the internal fractal dimension d of the dynamical graphs. We use the most naive 
definition of distance (the fewest links between two points ) so we are considering the "mathematical 
geometry" rather than the "physical geometry" in the terminology of [E5| R. We measured the internal 
fractal dimension at Pc for all the simulations, obtaining the values in the last column of Table 1 and then 
extrapolated to A^ = oo. There is however some ambiguity in this extrapolation: for Potts models P, p[ 
we simply extrapolate as N~^l'^'^ where vd is obtained from Binder's cumulant. Here we cannot do this 
so must assume a value for vd^ do the extrapolation and check that the value of d obtained is consistent. 
As c = 1 for both the (5 = 4 Potts and XY models, we use the value oi vd — 2 reported for the former 
in PJ. Then extrapolating d versus N^^l"^ yields intercept 2.66(3), very close to d = 2.7 which is usually 
obtained P, pi, p5|, p^. However with v — 1/2, which is what we have assumed throughout this analysis, 
vd = 1.33 rather than 2. Therefore we repeat the fit with vd = 1.33 and obtain d = 2.60(3). This is now 
consistent and we report our best estimate of the internal fractal dimension for the XY model coupled to 
2d quantum gravity as d = 2.6(1). 

4 Conclusions 

Our numerical results from simulations of the XY model on dynamical random phi-cubed graphs are 
consistent with the persistence of a KT type phase transition when the XY model is coupled to two- 
dimensional quantum gravity, in agreement with the theoretical predictions of [ [L7[ . We found that the 
critical point moved to lower temperature (from /3c = 2.3(1) to 2.4(1)), and the peak in the specific 
heat appeared to be growing more than on a fixed graph, but that the latter is probably due to modest 
statistics. The simulation on a single fixed graph generated in a pure two dimensional quantum gravity 
simulation, as for the fsing model |Q, gives results consistent with a regular lattice even though the 
status of the continuum limit in this case is uncertain. The universality in the graph properties we 
discovered in our earlier work on various Potts models with the same value of the central charge does not 
appear to carry over to the XY model, as it does not match the other c = 1 models we have simulated. 
A possible extension of this work would be to couple multiple XY models to quantum gravity and 
see what happens as c increases beyond 1. Similar investigations using multiple Potts models found no 
dramatic change B, but as the XY model does not appear to share the universal graph properties we 
found for the Potts models its behavior for c > 1 may be different. It would also be interesting to use the 
methods outlined in this paper to investigate the Villain form of the XY model, and higher dimensional 
0{N) spin models, coupled to quantum gravity. 
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^It now appears that the two are, in fact, identical ]27|. 
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Figure Captions 

Fig. 1. C for fixed simulations; points are from fluctuations in energy, line is from numerical differenti- 
ation of energy for N = 1000 only. 

Fig. 2. [ln{x) — IniC)] / ln{N)] = (2 — ri)/d versus (3 for fixed simulations. 

Fig. 3. [ln{x) — ln{C)] / ln{N)] = (2 — ?y)/<i versus /3 for dynamical simulations. 

Fig. 4. C for N = 100, 500, 1000 dynamical simulations; points are from fluctuations in energy, lines 
from numerical differentiation of energy. 

Fig. 5. AF, AL and PR3 for N ~ 2000 simulation; the y-scale applies to AF only, AL and PR3 have 
been scaled appropriately to flt on plot; f3c is indicated by vertical line. 

Fig. 6. PR3 for XY model along with Q = 2,3, 4, 10 Potts models (from M) for comparison; lines mark 
the Pc's. 



